function svd_applications()
% SVD_APPLICATIONS 奇异值分解的应用
% 
% 展示SVD在各种实际问题中的应用

fprintf('=== 奇异值分解应用 ===\n\n');

%% 1. 最小二乘问题求解
fprintf('1. 最小二乘问题求解\n');
fprintf('目标：使用SVD求解过定和欠定系统\n\n');

% 过定系统
m1 = 8; n1 = 4;
rng(42);
A1 = randn(m1, n1);
x1_true = [1; -2; 3; 0.5];
noise1 = 0.1 * randn(m1, 1);
b1 = A1 * x1_true + noise1;

fprintf('过定系统: %d×%d 矩阵\n', m1, n1);
fprintf('真实解: [%.2f, %.2f, %.2f, %.2f]''\n', x1_true);

% SVD分解
[U1, S1, V1] = svd(A1, 'econ');
singular_values1 = diag(S1);
fprintf('奇异值: [%.4f, %.4f, %.4f, %.4f]\n', singular_values1);
fprintf('条件数: %.2e\n', singular_values1(1) / singular_values1(end));

% SVD求解最小二乘
x1_svd = V1 * (S1 \ (U1' * b1));
fprintf('SVD求解结果: [%.4f, %.4f, %.4f, %.4f]''\n', x1_svd);

% 与其他方法比较
x1_normal = (A1' * A1) \ (A1' * b1);  % 正规方程
[Q1, R1] = qr(A1, 0);
x1_qr = R1 \ (Q1' * b1);              % QR分解

fprintf('正规方程结果: [%.4f, %.4f, %.4f, %.4f]''\n', x1_normal);
fprintf('QR分解结果:   [%.4f, %.4f, %.4f, %.4f]''\n', x1_qr);

% 残差比较
residual_svd = norm(A1 * x1_svd - b1);
residual_normal = norm(A1 * x1_normal - b1);
residual_qr = norm(A1 * x1_qr - b1);

fprintf('残差比较: SVD = %.2e, 正规方程 = %.2e, QR = %.2e\n', ...
        residual_svd, residual_normal, residual_qr);

%% 2. 欠定系统的最小范数解
fprintf('\n2. 欠定系统的最小范数解\n');
fprintf('目标：求解欠定系统的最小范数解\n\n');

% 欠定系统
m2 = 3; n2 = 5;
A2 = randn(m2, n2);
b2 = randn(m2, 1);

fprintf('欠定系统: %d×%d 矩阵\n', m2, n2);
fprintf('系统有无穷多解，寻找最小范数解\n');

% SVD分解
[U2, S2, V2] = svd(A2);
singular_values2 = diag(S2);
rank_A2 = sum(singular_values2 > 1e-12);

fprintf('矩阵的秩: %d\n', rank_A2);
fprintf('奇异值: [%.4f, %.4f, %.4f]\n', singular_values2(1:3));

% 最小范数解：x = V * S^+ * U' * b
S2_pinv = zeros(size(S2'));
for i = 1:rank_A2
    S2_pinv(i, i) = 1 / singular_values2(i);
end

x2_min_norm = V2 * S2_pinv * U2' * b2;
fprintf('最小范数解的范数: %.4f\n', norm(x2_min_norm));

% 验证解的正确性
residual2 = norm(A2 * x2_min_norm - b2);
fprintf('残差: %.2e\n', residual2);

% 与MATLAB内置函数比较
x2_pinv = pinv(A2) * b2;
fprintf('MATLAB pinv结果的范数: %.4f\n', norm(x2_pinv));
fprintf('两种方法的差异: %.2e\n', norm(x2_min_norm - x2_pinv));

%% 3. 矩阵的低秩近似
fprintf('\n3. 矩阵的低秩近似\n');
fprintf('目标：使用SVD进行矩阵的低秩近似\n\n');

% 构造低秩矩阵加噪声
rank_true = 3;
m3 = 10; n3 = 8;
U3_true = orth(randn(m3, rank_true));
V3_true = orth(randn(n3, rank_true));
S3_true = diag([10, 5, 2]);
A3_clean = U3_true * S3_true * V3_true';

% 添加噪声
noise_level = 0.5;
A3_noisy = A3_clean + noise_level * randn(m3, n3);

fprintf('原始矩阵: %d×%d，真实秩 = %d\n', m3, n3, rank_true);
fprintf('噪声水平: %.1f\n', noise_level);

% SVD分解
[U3, S3, V3] = svd(A3_noisy);
singular_values3 = diag(S3);

fprintf('奇异值: [%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f]\n', ...
        singular_values3);

% 不同秩的近似
ranks_to_test = [1, 2, 3, 4, 5];
fprintf('\n低秩近似分析:\n');
fprintf('秩    近似误差    压缩比    相对误差\n');
fprintf('--    --------    ------    --------\n');

for r = ranks_to_test
    % r秩近似
    A3_approx = U3(:, 1:r) * S3(1:r, 1:r) * V3(:, 1:r)';
    
    % 近似误差
    approx_error = norm(A3_noisy - A3_approx, 'fro');
    relative_error = approx_error / norm(A3_noisy, 'fro');
    
    % 压缩比
    original_elements = m3 * n3;
    compressed_elements = r * (m3 + n3 + 1);
    compression_ratio = compressed_elements / original_elements;
    
    fprintf('%2d    %8.4f    %6.3f    %8.4f\n', ...
            r, approx_error, compression_ratio, relative_error);
end

% 最优秩的选择（基于奇异值的跳跃）
singular_ratios = singular_values3(1:end-1) ./ singular_values3(2:end);
[~, optimal_rank] = max(singular_ratios);
fprintf('建议的最优秩: %d (基于奇异值比值)\n', optimal_rank);

%% 4. 主成分分析 (PCA)
fprintf('\n4. 主成分分析 (PCA)\n');
fprintf('目标：使用SVD进行数据降维\n\n');

% 生成相关数据
n_samples = 100;
n_features = 5;
rng(123);

% 构造具有内在低维结构的数据
true_components = 3;
W_true = randn(n_features, true_components);
H_true = randn(true_components, n_samples);
X4_clean = W_true * H_true;
X4_noisy = X4_clean + 0.5 * randn(n_features, n_samples);

fprintf('数据维度: %d×%d\n', size(X4_noisy));
fprintf('真实内在维度: %d\n', true_components);

% 数据中心化
X4_mean = mean(X4_noisy, 2);
X4_centered = X4_noisy - X4_mean;

% SVD分解进行PCA
[U4, S4, V4] = svd(X4_centered, 'econ');
singular_values4 = diag(S4);

% 主成分的方差解释
total_variance = sum(singular_values4.^2);
explained_variance_ratio = singular_values4.^2 / total_variance;
cumulative_variance = cumsum(explained_variance_ratio);

fprintf('主成分分析结果:\n');
fprintf('成分    奇异值    方差解释比    累积解释比\n');
fprintf('----    ------    ----------    ----------\n');
for i = 1:min(5, length(singular_values4))
    fprintf('%4d    %6.2f    %10.3f    %10.3f\n', ...
            i, singular_values4(i), explained_variance_ratio(i), cumulative_variance(i));
end

% 选择主成分数量（保留95%方差）
n_components = find(cumulative_variance >= 0.95, 1);
fprintf('保留95%%方差需要 %d 个主成分\n', n_components);

% 降维和重构
X4_reduced = S4(1:n_components, 1:n_components) * V4(:, 1:n_components)';
X4_reconstructed = U4(:, 1:n_components) * X4_reduced + X4_mean;

% 重构误差
reconstruction_error = norm(X4_noisy - X4_reconstructed, 'fro') / norm(X4_noisy, 'fro');
fprintf('重构相对误差: %.4f\n', reconstruction_error);

%% 5. 图像压缩
fprintf('\n5. 图像压缩\n');
fprintf('目标：使用SVD压缩图像数据\n\n');

% 创建简单的"图像"矩阵
img_size = 50;
[X_img, Y_img] = meshgrid(1:img_size, 1:img_size);
I5 = sin(2*pi*X_img/img_size) .* cos(2*pi*Y_img/img_size) + ...
     0.5 * sin(4*pi*X_img/img_size) + 0.3 * randn(img_size, img_size);

fprintf('图像尺寸: %d×%d\n', size(I5));

% SVD分解
[U5, S5, V5] = svd(I5);
singular_values5 = diag(S5);

fprintf('前10个奇异值: [%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f]\n', ...
        singular_values5(1:10));

% 不同压缩级别
compression_ranks = [1, 5, 10, 20, 30];
fprintf('\n图像压缩分析:\n');
fprintf('秩    压缩比    PSNR(dB)    相对误差\n');
fprintf('--    ------    --------    --------\n');

for r = compression_ranks
    % 压缩重构
    I5_compressed = U5(:, 1:r) * S5(1:r, 1:r) * V5(:, 1:r)';
    
    % 压缩比
    original_size = numel(I5);
    compressed_size = r * (size(I5, 1) + size(I5, 2) + 1);
    compression_ratio = compressed_size / original_size;
    
    % 图像质量评估
    mse = mean((I5(:) - I5_compressed(:)).^2);
    psnr = 10 * log10(max(I5(:))^2 / mse);
    relative_error = norm(I5 - I5_compressed, 'fro') / norm(I5, 'fro');
    
    fprintf('%2d    %6.3f    %8.2f    %8.4f\n', ...
            r, compression_ratio, psnr, relative_error);
end

%% 6. 推荐系统
fprintf('\n6. 推荐系统\n');
fprintf('目标：使用SVD进行协同过滤\n\n');

% 模拟用户-物品评分矩阵（稀疏）
n_users = 20;
n_items = 15;
sparsity = 0.3;  % 30%的评分已知

% 生成潜在因子模型
n_factors = 4;
U6_true = 5 * rand(n_users, n_factors);  % 用户特征
V6_true = 5 * rand(n_items, n_factors);  % 物品特征
R6_true = U6_true * V6_true';            % 真实评分矩阵

% 创建稀疏观测
observed_mask = rand(n_users, n_items) < sparsity;
R6_observed = R6_true;
R6_observed(~observed_mask) = 0;  % 未观测的设为0

fprintf('评分矩阵: %d用户 × %d物品\n', n_users, n_items);
fprintf('观测密度: %.1f%%\n', 100 * sparsity);
fprintf('观测评分数: %d / %d\n', sum(observed_mask(:)), numel(observed_mask));

% 简单的SVD方法（将0当作真实评分，实际中需要更复杂的方法）
[U6, S6, V6] = svd(R6_observed);
singular_values6 = diag(S6);

fprintf('前5个奇异值: [%.2f, %.2f, %.2f, %.2f, %.2f]\n', ...
        singular_values6(1:5));

% 不同因子数的重构
factors_to_test = [2, 4, 6, 8];
fprintf('\n协同过滤分析:\n');
fprintf('因子数    观测误差    预测误差    总体误差\n');
fprintf('------    --------    --------    --------\n');

for k = factors_to_test
    % k因子近似
    R6_approx = U6(:, 1:k) * S6(1:k, 1:k) * V6(:, 1:k)';
    
    % 观测位置的误差
    observed_error = norm((R6_true - R6_approx) .* observed_mask, 'fro') / ...
                     norm(R6_true .* observed_mask, 'fro');
    
    % 未观测位置的误差（预测误差）
    unobserved_mask = ~observed_mask;
    prediction_error = norm((R6_true - R6_approx) .* unobserved_mask, 'fro') / ...
                       norm(R6_true .* unobserved_mask, 'fro');
    
    % 总体误差
    total_error = norm(R6_true - R6_approx, 'fro') / norm(R6_true, 'fro');
    
    fprintf('%6d    %8.4f    %8.4f    %8.4f\n', ...
            k, observed_error, prediction_error, total_error);
end

%% 7. 信号处理中的去噪
fprintf('\n7. 信号处理中的去噪\n');
fprintf('目标：使用SVD进行信号去噪\n\n');

% 构造信号矩阵（Hankel矩阵）
t = 0:0.01:2*pi;
signal_clean = sin(2*t) + 0.5*sin(5*t) + 0.3*cos(3*t);
signal_noisy = signal_clean + 0.3 * randn(size(t));

fprintf('信号长度: %d 点\n', length(t));
fprintf('信噪比: %.2f dB\n', 20*log10(std(signal_clean)/std(signal_noisy - signal_clean)));

% 构造Hankel矩阵
window_size = 100;
n_windows = length(signal_noisy) - window_size + 1;
H7 = zeros(window_size, n_windows);

for i = 1:n_windows
    H7(:, i) = signal_noisy(i:i+window_size-1);
end

fprintf('Hankel矩阵尺寸: %d×%d\n', size(H7));

% SVD分解
[U7, S7, V7] = svd(H7);
singular_values7 = diag(S7);

fprintf('前10个奇异值: [%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f]\n', ...
        singular_values7(1:10));

% 选择主要成分进行去噪
n_signal_components = 6;  % 基于信号的频率成分
H7_denoised = U7(:, 1:n_signal_components) * S7(1:n_signal_components, 1:n_signal_components) * ...
              V7(:, 1:n_signal_components)';

% 从Hankel矩阵重构信号（取对角线平均）
signal_denoised = zeros(size(signal_noisy));
for i = 1:length(signal_denoised)
    if i <= window_size
        signal_denoised(i) = mean(diag(H7_denoised, i-1));
    else
        signal_denoised(i) = mean(diag(H7_denoised, i-window_size));
    end
end

% 简化重构：只取第一行
signal_denoised = H7_denoised(1, :);
if length(signal_denoised) < length(signal_noisy)
    signal_denoised = [signal_denoised, signal_noisy(length(signal_denoised)+1:end)];
end

% 去噪效果评估
snr_original = 20*log10(norm(signal_clean)/norm(signal_noisy - signal_clean));
snr_denoised = 20*log10(norm(signal_clean(1:length(signal_denoised)))/...
                        norm(signal_clean(1:length(signal_denoised)) - signal_denoised));

fprintf('去噪效果:\n');
fprintf('原始信噪比: %.2f dB\n', snr_original);
fprintf('去噪后信噪比: %.2f dB\n', snr_denoised);
fprintf('改善: %.2f dB\n', snr_denoised - snr_original);

%% 8. 数值秩的确定
fprintf('\n8. 数值秩的确定\n');
fprintf('目标：使用SVD确定矩阵的数值秩\n\n');

% 构造不同数值秩的矩阵
test_matrices = {
    randn(10, 8),                           % 满秩矩阵
    randn(10, 5) * randn(5, 8),            % 秩5矩阵
    randn(10, 3) * randn(3, 8) + 1e-10 * randn(10, 8)  % 秩3矩阵+噪声
};

matrix_names = {'满秩矩阵', '秩5矩阵', '秩3矩阵+噪声'};

fprintf('数值秩分析:\n');
for i = 1:length(test_matrices)
    A_test = test_matrices{i};
    [~, S_test, ~] = svd(A_test);
    sv_test = diag(S_test);
    
    fprintf('\n%s (%d×%d):\n', matrix_names{i}, size(A_test));
    fprintf('奇异值: [%.2e, %.2e, %.2e, %.2e, %.2e, ...]\n', sv_test(1:5));
    
    % 不同阈值下的数值秩
    tolerances = [1e-12, 1e-10, 1e-8, 1e-6, 1e-4];
    fprintf('阈值        数值秩\n');
    fprintf('--------    ------\n');
    
    for tol = tolerances
        numerical_rank = sum(sv_test > tol);
        fprintf('%.0e    %6d\n', tol, numerical_rank);
    end
    
    % MATLAB的rank函数结果
    matlab_rank = rank(A_test);
    fprintf('MATLAB rank: %d\n', matlab_rank);
end

%% 总结
fprintf('\n=== SVD应用总结 ===\n');
fprintf('1. 最小二乘求解：处理过定和欠定系统\n');
fprintf('2. 低秩近似：数据压缩和去噪\n');
fprintf('3. 主成分分析：数据降维和特征提取\n');
fprintf('4. 图像压缩：有损压缩的数学基础\n');
fprintf('5. 推荐系统：协同过滤和矩阵补全\n');
fprintf('6. 信号去噪：基于子空间的方法\n');
fprintf('7. 数值秩确定：处理数值误差\n\n');

fprintf('SVD的优势：\n');
fprintf('• 数值稳定性好\n');
fprintf('• 适用于任意矩阵\n');
fprintf('• 提供最优低秩近似\n');
fprintf('• 几何意义清晰\n\n');

fprintf('实际应用考虑：\n');
fprintf('• 计算复杂度较高\n');
fprintf('• 需要选择合适的秩\n');
fprintf('• 处理稀疏矩阵的特殊方法\n');
fprintf('• 大规模问题的近似算法\n\n');

end